Two-loop QED corrections with closed fermion loops 



00 . 
O ■ 

o ■ 

<n : 

5-h ■ 
Oh' 

<: 

m '. 

(N ■ 



: 

I , 



: 

O 1 

• i-H ' 

C/5 ■ 
>Y 

^ : 

Ph. 



> 

in 
o 

o 

oo 
o 



V. A. Yerokhin, 1 P. Indelicato, 2 and V. M. Shabaev 3 

J Center for Advanced Studies, St. Petersburg State Polytechmcal University, 
Polytekhnicheskaya 29, St. Petersburg 195251, Russia 
2 Laboratoire Kastler-Brossel, Ecole Normale Superieure, 
CNRS et Universite P. et M. Curie, Case 74, 4 pi. Jussieu, F-75252, France 
s Department of Physics, St. Petersburg State University, Oulianovskaya 1, St. Petersburg 198504, Russia 

We report a calculation of all two-loop QED corrections with closed fermion loops for the n = 1 
and n = 2 states of H-like ions and for a wide range of the nuclear charge numbers Z = 1 — 100. The 
calculation is performed to all orders in the binding-strength parameter Za, with the exception that 
in a few cases the free-loop approximation is employed in the treatment of the fermion loops. Detailed 
comparison is made with previous Za-expansion calculations and the higher-order remainder term 
to order a 2 (Za) 6 is identified. 
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Introduction 

Highly charged ions are often considered as the ideal 
testing ground for investigating the strong-field regime of 
bound-state quantum electrodynamics (QED). They fea- 
ture a strong static Coulomb field of the nucleus and have 
a simple electronic structure, which can be accurately 
determined in ab initio theoretical calculations. In this 
respect, the ultimate investigation object is the H-likc 
uranium, the heaviest naturally occurring element. Mea- 
surements of the Is Lamb shift in H-like uranium have 
progressed drastically during last decades [l|, , having 
achieved an accuracy of 4.6 eV [H, which corresponds 
to a fractional error of 1.7% with respect to the total 
QED contribution. Further advance anticipated in the 
future will make such experiments sensitive to the two- 
loop QED effects. 

Even higher precision is achieved for heavy Li-like 
ions. Measurements of the 2p 1 / 2 ,3/2-2s transition ener- 
gies [3, H, H 0, HI have lately reached a fractional accu- 
racy of 0.03% with respect to the total QED contribution. 
This corresponds to a 10% sensitivity of the experimental 
results to the two-loop QED effects. These measurements 
provide an excellent possibility for the identification of 
the two-loop Lamb shift and for testing the bound-state 
QED up to second order in a in the strong-field regime (a 
is the fine-structure constant). Theoretical description of 
Li-like ions is more complicated than that of H-like ions, 
which is due to the presence of additional electrons. For 
heavy systems, the electron-electron interaction is weak 
(as compared to the electron-nucleus interaction) and can 
be successfully treated by perturbation theory. Ab initio 
calculational results for the effects of the electron corre- 
lation and the screening of one-loop QED corrections are 
already available for Li- like ions [3, [lCj ; their accuracy is 
sufficient for identification of the two-loop QED effects. 

The detailed theoretical understanding of the two-loop 
QED effects is also necessary for the interpretation of 
high-precision experimental data in the low-Z region. 
The most prominent example here is hydrogen. Its spec- 



troscopy can nowadays be realized with a relative ac- 
curacy on the level of 10" 14 p], El- The theoretical 
understanding of the Is and 2s Lamb shift in hydrogen 
is still limited, to a large extent, by the two-loop QED 
effects [13l | . 

The subject of our present interest is the set of the two- 
loop QED corrections (also referred to as the two-loop 
Lamb shift), graphically represented in Fig. [TJ These 
corrections have been intensively investigated within the 
perturbative expansion in the binding-strength parame- 
ter Za fLU [T5L fl6l. [17(1 (Z is the nuclear charge number). 
The results of these studies, however, do not provide re- 
liable information about the magnitude of the two-loop 
Lamb shift in heavy ions, where the parameter Za ap- 
proaches unity. A non-perturbative (in Za) evaluation 
of the whole set of the two-loop diagrams is needed for 
the interpretation of experimental results available in the 
middle- and high-Z region. 

Numerical all-order calculations of the two-loop cor- 
rections started in late 1980s [H, QJ, HO, HH E3, i| and 
dealt with the diagrams with the closed fermion loops 
[Fig. [ljd)-(k)]. An evaluation of the three remaining di- 
agrams [Fig. Dla)-(c)], referred to as the two-loop self- 
energy diagrams, turned out to require considerable ef- 
forts. It was accom plis hed in a series of investigations 
0, [H, [H, H3, HI, IM H3 for the nuclear charge 
numbers Z > 10 for the ground state and Z > 60 for the 
n = 2 states. 

The goal of this work is to perform a detailed investi- 
gation of all two-loop diagrams with the closed fermion 
loops [Fig. Qld)-(k)], extending previous evaluations to 
cover the whole region of the nuclear charge numbers 
Z = 1 — 100 and all n = 1 and n = 2 states. The first 
results of our calculation were presented in Ref . [3l| . At 
present, our intention is to achieve high numerical ac- 
curacy in the \ow-Z region. This will allow us to make 
a detailed comparison with the perturbative calculations 
and to isolate the non-perturbative remainder to order 
a 2 (Za) 6 , which is of experimental interest for hydrogen. 
Our calculation will be performed to all orders in the 
binding-strength parameter Za, but an approximation 
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will be made in the treatment of the fermion loops in the 
diagrams in Fig. HJh)-(k). This approximation, referred 
to as the free-loop one, consists in keeping the leading 
nonvanishing term of the expansion of the fermion loop 
in terms of the binding potential. In the one-loop case, it 
corresponds to the Uchling potential and yields the dom- 
inant contribution even for high-Z ions. To perform a 
non-perturbative calculation of the diagrams in Fig.fljh)- 
(k) without the free-loop approximation is a difficult and 
a so far unsolved problem. 

The paper is organized as follows. In Section HI we 
describe our approach to the evaluation of the one-loop 
self-energy and vacuum-polarization corrections, which is 
used as a basis for treatment of the two-loop corrections. 
In Section m we discuss the evaluation of individual two- 
loop corrections with closed fermion loops, present our 
numerical results, and compare them with the data ob- 
tained within the perturbative Za-expansion approach. 
In Section [lll| we summarize our results. 

The relativistic units (h = c = 1) are used in this 
paper. 



I. ONE-LOOP QED CORRECTIONS 

We start with presenting some basic formulas for the 
first-order self-energy and vacuum-polarization correc- 
tions that will be needed below in addressing the two- 
loop corrections. 




(d) (e) (f) 




(g) (h) (i) (k) 



FIG. 1: Two-loop one-electron QED corrections. Gauge- 
invariant subsets are referred to as SESE (a-c), SEVP (d-f), 
VPVP (g-i), and S(VP)E (k). 



A. Self-energy 

The one-loop self-energy (SE) correction to the Lamb 
shift can be represented as a matrix element of the renor- 
malized SE operator E R , 

AE SE = (a| 7 £ R (£ a )|a) , (1) 

where E^ = E — 8m, 8m is the one-loop mass coun- 
terterm, e a is the energy of the reference state, E is the 
unrcnormalized SE operator, 

/oo 
-oo 

xG(£-w,Si,x 2 )q f D^(w,a;i2) ,(2) 

G is the Dirac Coulomb Green function G(e) = [e— 7i{l — 
i0)]~ , H is the Dirac Coulomb Hamiltonian, D^ v is the 
photon propagator, a} 1 = (1, a), and X12 = aci — X2- 
It is assumed that the above expressions are regularized 
in a covariant way and that the limit that removes the 
regularization is properly approached. 

In order to facilitate the numerical evaluation of the 
above expressions, one needs to explicitly eliminate di- 
vergences from the SE operator. A popular method of 
doing this [32j is based on the expansion of the SE oper- 
ator in terms the binding field, 

E = E(°)+E( 1 >+E( 2 +\ (3) 

where the superscript (i), i — 0,1 indicates the total 
number of interactions with the binding Coulomb field 
and the index (2+) labels the term generated by > 2 
such interactions. Only the first two terms of the ex- 
pansion are divergent; all divergences can be shown to 
vanish in the difference E — 8m. The divergent parts are 
regularized by working in an extended number of dimen- 
sions and evaluated in momentum space, see Ref. [3 31] for 
details. We mention that the first-order expansion term 
E 1 - 1 ) is usually represented as a product of the time com- 
ponent of the standard vertex operator T a (p\,p2) and 
the Coulomb potential Vc, 

= r°( Pl ,p 2 ) V c (\pi -pa|) . (4) 

The energy shifts induced by (the final parts of) the 
three terms in the right-hand-side of Eq. ((3J are referred 
to as the zero-, one-, and many-potential terms, respec- 
tively: 

AE SE = AE^> + A.EgE + A4 2 E +) , (5) 

where 

a4e = / j^s 4%) MP) , (6) 

(7) 
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Ae se )= J dx 1 dx 2 ^ a {x l )^ 2+ \e a ,x ll x 2 )^ a {x 2 ) 1 

(8) 

where q = Pi - p 2 , tp a = *Pll°, V C (q) = -AnZa/\q\ 2 
is the Coulomb potential in momentum space, E v (p) 
is the finite part of the subtracted free SE operator 
j](o) (p^ — firn, r^(pi , p 2 ) is the finite part of the vertex op- 
erator r /i (pi,p2); and p, p\, and p 2 are four- vectors with 
the time component fixed by p a = p\ = p 2 = e a . The op- 
erator £*- 2+ ' is obtained from Eq. ((2]) by the substitution 
G — > G^ 2+ \ where the index (2+) denotes, as usual, two 
or more interactions with the binding Coulomb field. 

In the extended number of dimensions, D = 4 — 2e, the 
free SE operator £(°) is given by 



£ (0) O) = -Aniap 2 



d D k 1 J a (#- ^ + m)7 CT 



(2ir) D k 2 (p - k) 2 - m 2 



(9) 

where j> = p v "i v and p is the auxiliary mass parameter 
introduced in order to keep the proper dimension of the 
interaction term in the Lagrangian. The momentum in- 
tegration in Eq. © is simple; it is described, e.g., in 
Appendix A of Ref. [28| . Omitting terms of order e and 
higher, one obtains 



£(°)(p)-£m = 



aC e 



a 
An 



2m 1 



1-P 



2/Q 

1-p 

P 

1-P 



In p 
lnp 



(10) 



where p = (m 2 — p 2 )/m 2 and G c = T(l + 
e)(47r) e (p 2 /m 2 )\ The operator in Eq. © is the 
finite part of the right-hand-side of Eq. ([1^) when e ap- 
proaches zero. 

The free vertex operator T a is 



F a (pi,p 2 ) = — Airiap 

X7 



d D fc 1 *4 a - U 



(2tt) d fc 2 Ox - fc) 2 - m 2 
i> 2 - $ + m 



(p 2 — fc) 2 — m 2 
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(11) 



The momentum integration is performed by usin g th e 
standard technique (see, e.g., Appendix A of Ref. [28j| ) . 
Omitting terms of order e and higher, one obtains 



aC e a n 



2 7 



n 



12) 



The operator T^, in Eq. ([6]) is the finite part of the right- 
hand-side of Eq. (|12[) when e approaches zero. 

Out of the three terms in Eq. ([5]), the third one 
[A_Eg E ] is the most difficult to evaluate numerically. To 
a large extent, this is due to the partial-wave expansion 
that inevitably appears in the evaluation of the function 
G^ 2+ \ The convergence of this expansion is rather good 
for the ground state of high-Z ions (provided that ra- 
dial integrations are carried out first), but quickly wors- 
ens when Z is decreased and(or) the principal quantum 
number n is increased. An approach to overcome this 
difficulty was suggested in Ref. [34] . It is based on the 
separation of the > 2-potential Green function G^ 2-1- ) into 
two parts, 

G ( - 2+ \e,x 1) x 2 ) =G (2+ \e, Xl ,x 2 ) 

+ [G^(e, Xl ,x 2 )-G^(s,x u x 2 )] , 

(13) 

where G a 2+ ^ yields an approximation to G^ 2+ ) in the re- 
gion where X\ ss x 2 and the energy argument is far away 
from a pole. This separation is related to the one used 
in Ref. [35| for the identification of ultraviolet divergent 
parts of the SE operator in coordinate space. The func- 

(2+) 

tion G a is given by 

Gi 2+ ) (e, x u x 2 ) = G<® (e + fi, x u x 2 ) - G<°) (e, x x ,x 2 ) 
d 



(14) 



E=e 



where G^ is the free Dirac Green function and fl = 

(2+) 

2Za/(xi + x 2 ). The function G a has two important 
properties: (i) it can be easily expressed both in the 
closed form and in the partial-wave expansion form and 
(ii) it depends on angular variables through X\ 2 = X\—x 2 
only. These features ensure that the numerical evaluation 

(2+) 

of expressions with G a is rather simple. The main rea- 

(2+) 

son to introduce the auxiliary function Ga is that the 
separation (|13p improves the convergence of the partial- 
wave expansion. It was demonstrated in Ref. [34j that 
the substitution of the difference 

G (2+) _ q(2+) intQ (thc 

high-energy part of) Eq. ([5]) improves the convergence of 
the partial-wave expansion drastically and allows one to 
obtain reasonably accurate results for the SE correction 
at Z = 1 with employing just about 30 partial waves. 
This approach was extensively used in the present work 
for the evaluation of the two- loop corrections. Its usage 
allowed us to obtain accurate numerical results for the 
whole region of the nuclear charge numbers Z = 1 — 100 
in calculations of the diagrams in Fig. [ijd)-(f) and (k). 



B. Vacuum-polarization 



where N a (xb) = 7 CT (^i — xf> + m)-f a (^) 2 — xfi + m)^, 
b = ypi + {l — y)p2, and D = xb 2 + m 2 - yp\ - (1 - y)p\. 



The one-loop vacuum-polarization (VP) potential con- 
sists of two parts, which are commonly referred to as the 
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Uehling and the Wichmann-Kroll ones, see the review 
HI] for details, 

U VP (r) = C/ Uch (r) + U WK (r) . (15) 

The expression for the Uehling potential is given by 

2a 2 Z 



U\jch(r) = - 



3mr 



x [K (2m\r - r'\) - K (2m\r + r'\)} , (16) 



where 

K (x) = 



dt( 



1 

2fi 



(17) 



and the nuclear-charge density p(r) is spherically sym- 
metric and normalized by the condition J dr p(r) = 1. 
The Wichmann-Kroll potential is conveniently expressed 
in terms of the VP charge density py P (r), 



U WK (r) 



/ dr'r' 2 -p^\r') 
Jo r > 



(18) 



where r> = 
written as 



^(r) = ^Re 

7T 



(3+) 

where G K denotes the radial part of the electron propa- 
gator that contains three and more interactions with the 
binding Coulomb field and k is the relativistic angular 
quantum number. It was shown [3?l |38| that no spu- 
rious terms arise in a numerical evaluation of Eq. (fT5|) 
if the expansion over k is terminated by a finite cutoff 
parameter. 

According to the Furry theorem, all parts of the elec- 
tron propagator that contain an even number of interac- 
tions with the binding field yield a vanishing contribution 
to the VP charge density. Thus, in numerical evalua- 
tions, the function Gl 3+ ^ in Eq. (TT§|) can be substituted 



max(r, r'). The VP charge density can be 



E 



(19) 



(2+) 

by G\ , which is conveniently expressed as 
Gl 2 +) (w, x, y) = / dz 7? Gf (w, x, z) V(z) 



x [G K (u,z,y)-G<£\cj,z,y)}, (20) 

where G K is the radial part of the bound-electron prop- 
agator, G« is the radial part of the free electron propa- 
gator, and V(z) is the binding potential. 

Formulas (|16p -(|20 p were employed for the numerical 
evaluation of the one-loop VP potential in this work us- 
ing the numerical procedure developed in Refs. [H, H^]. 
The Uehling potential was calculated with the Fermi 
model of the nuclear-charge distribution, whereas for 
the Wichmann-Kroll potential, the spherical shell model 



[p(r) oc 8(r — R)] was employed. The summation over k 
in Eq. (JTUJ) was extended up to |K max | = 10. 

We mention that in the case when high numerical ac- 
curacy is not needed and the finite nuclear size correction 
may be disregarded, the Wichmann-Kroll potential can 
be conveniently evaluated by employing analytical ap- 
proximation formulas reported in Ref. [401 ] . 



II. TWO-LOOP QED CORRECTIONS 

The two-loop contributions to the energy shift are con- 
veniently represented in terms of the dimensionless func- 
tion F(Za) defined by 



/a\2 (Za) 4 
A.E — m - i — 
\7r/ rr 



. F(Za), 

i 

where n is the principal quantum number. 



(21) 



A. Self-energy in the Coulomb potential modified 
by the vacuum polarization 

We start our consideration of the two-loop corrections 
with the set of the three diagrams in Fig. [ljd)-(f). This 
set is gauge invariant and can be regarded as the one- 
loop SE correction in the combined field of the nuclear 
Coulomb and the VP potential. The corresponding en- 
ergy shift will be referred to as the SEVP correction in 
the following. 

Formal expressions for the SEVP correction can be ob- 
tained by considering the first-order perturbation of the 
one-loop SE correction by the VP potential (fTS"]) . Pertur- 
bations of the reference-state wave function, the binding 
energy, and the electron propagator give rise to the ir- 
reducible, the reducible, and the vertex contributions, 
respectively. The irreducible part is 



(a| 7 E i j(e a )|5a) + (Sa\^ a (e a )\a) , (22) 



where E# is the renormalized SE operator, and \Sa) is 
the first-order perturbation of the reference-state wave 
function \a) by Uyp. The reducible part is given by 



A£?evp = (a\U VP \a) (a| 7o ^E fl ( £ )|a) 
and the vertex part is 



(23) 



^^SEVP 



2ia 



nin 2 



(ni\Uyp\n 2 ) {an 2 \a D 1 - 1 " (uj)\nia) 
(e a -u - e ni )(e a — uj - e„ 2 ) 



(24) 



where the summation over ni^ goes over the Dirac spec- 
trum, and the virtual-state energies are assumed to have 
a small imaginary addition, e n — > e„(l — i0). 
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The problem of calculating the SE correction in the 
presence of the perturbing potential has been extensively 
studied in the literature, see, e.g., Ref. [HI, S3, EE S3 • In 
this work, we employ the general scheme developed in our 
previous study [45J] for the case of the SE correction to 
the hyperfine splitting. Several modifications were intro- 
duced into the scheme, among them the inclusion of the 
finite nuclear size effect. This modification was essential, 
firstly, because the effect is significantly enhanced by the 
singular behavior of the VP potential and, secondly, be- 
cause the extended nuclear charge distribution removes 
the logarithmic singularity of the point-nucleus Uchling 
potential, which simplifies numerical integrations consid- 
erably Even with the extended nuclear size, the usage of 
extremely fine grids was required in the nuclear region, 
in order to achieve a high controllable accuracy in radial 
integrations. 

In actual calculations, the Fermi model of the nuclear- 
charge distribution was employed for the evaluation of 
the reference-state wave function and the Uehling po- 
tential, whereas the electron propagator(s) inside the SE 
loop and the Wichmann-Kroll potential were calculated 
with the spherical shell model [p(r) cx 5{r — R)]. For 
systems with Z < 10, we neglected the nuclear-size de- 
pendence in the electron propagators. The values of the 
root-mean-square (rms) radii of the nuclear-charge dis- 
tribution were taken from Ref. [46] in most cases. For 

l li 

uranium, we used the value (r 2 ) = 5.8569 (33) fm ob- 
tained in the recent re-evaluation of experimental data 
[47l ]. In the case of fermium (Z = 100), there is no ex- 
perimental results available, so we used the interpolation 
formula from Ref. [48j and assigned a (conventional) un- 
certainty of 1% to the resulting value. The rms radii used 
in the present investigation are listed in second column 
of Table H 

The calculational results for the SEVP correction for 
the Is, 2s, 2p\/2, and 2p 3 / 2 states of H-like ions are pre- 
sented in Table [l] expressed in terms of the function 
F{Za) defined by Eq. (|2Tj) . Our results are in good 
agreement with the values reporte d p reviously for ura- 
nium, lead, and ytterbium in Ref. [201 ]. The uncertainty 
specified in the table includes the numerical error and 
the estimated errors due to the models of the nuclear 
charge distribution and due to uncertainties of the rms 
radii. The model dependence of the results was estimated 
by switching between the Fermi, the uniform, and the 
spherical-shell models in evaluations of the Uehling po- 
tential and the wave functions; the largest deviation was 
taken as the error due to the nuclear model. The er- 
ror due to the uncertainty of the nuclear radius was ob- 
tained by repeating the calculations with the rms radii 
varied within the error bars specified in the table. All 
three errors were added quadratically. In Table U and 
in the tables that follow, the omitted uncertainty means 
that the expected error is smaller than the last significant 
digit specified. 

It is instructive to compare our nonpertubative re- 
sults with the ones obtained within the Za-expansion 




10 20 30 40 50 60 70 80 90 100 

Nuclear charge number Z 



FIG. 2: Higher-order remainder G(Za) for the SEVP correc- 
tion. Errors due to the model of the nuclear-charge distribu- 
tion and the rms radii are not shown. 



approach. The Za expansion of the SEVP correction 
reads 

2 (Zrv) 5 



AEsevp =m ) - ^ |a 50 + {Za) In 2 [{Za) 2 ] a 62 



{Za) In [{Za)- 2 } a 61 + {Za) G(Za)} 



(25) 



where the function G{Za) is the higher-order remainder, 
G{0a) = ago . The results known for the coefficients of 

this expansion are: M M HE Hi: 



a 50 = 1.920576*1,0, 
4 x 

«62 = 77 Oi o , 

45 

061 = tjj (J + In 2 ) 



(26) 
(27) 



32/3 1 1 

45U + 4^-n +7 + * (n) - lnnKo 



n 2 - 1 



135 n 2 



(28) 



The higher-order remainder G{Za) was inferred from our 
numerical results for the SEVP correction and plotted on 
Fig. [2] For hydrogen, the results of our direct numerical 
evaluation are: G ls {a) = -13.2(4), G 2s {a) = -11.7(4), 
G 2pi/2 {a) = -0.034, and G 2p3/2 {a) = 0.015. For the 
normalized difference of the 2s and Is states and for the 
fine-structure difference, these values are consistent with 
the analytical results ago(2s) — O6o(ls) = 1.491199 and 
060(2^3/2) - aeo(2pi/ 2 ) = 1/20 [52|. 
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TABLE I: Energy shifts due the SEVP correction, in units of F(Za). Uncertainties specified include the estimated errors due 
to the nuclear charge distribution models and the values of the rms radius. 



z 


<r 2 ) 1/2 [fm] 


Is 


2s 


2pi/2 


2^3/2 


1 


0.879 (9) 


0.01474 (2) 


0.01459 (2) 


—0.000025 


—0.000022 


2 


1.676 (3) 


0.02993 (2) 


0.02945 (4) 


-0.000086 


-0.000076 


3 


2.43 (3) 


0.04519 (3) 


0.04430 (5) 


-0.000175 


-0.000152 


5 


2.41 (3) 


0.07546 (4) 


0.07355 (6) 


—0.000414 


-0.000352 


7 


2.558 (7) 


0.10508 (3) 


0.10208 (7) 


-0.000715 


-0.000597 


10 


3.005 (2) 


0.14804 (3) 


0.14351 (5) 


—0.001242 


-0.001008 


15 


3.189 (2) 


0.21618 (3) 


0.20961 (4) 


-0.002206 


-0.001700 


20 


3.476 (1) 


0.28099 (5) 


0.27394 (5) 


—0.003130 


-0.002265 


25 


3.706 (2) 


0.34394 (4) 


0.33841 (6) 


-0.003868 


-0.002577 


30 


3.929 (1) 


0.40649 (5) 


0.40497 (6) 


-0.004269 (1) 


-0.002529 (1) 


35 


4.163 (2) 


0.47012 (3) 


0.47570 (6) 


-0.004161 (1) 


-0.002025 (1) 


40 


4.270 (1) 


0.53642 (4) 


0.55290 (5) 


-0.003332 (1) 


-0.000974 (4) 


45 


4.494 (2) 


0.60687 (3) 


0.63898 (5) 


-0.001500 (1) 


0.000711 (2) 


50 


4.654 (1) 


0.68338 (3) 


0.73709 (5) 


0.001724 (2) 


0.003116 (6) 


55 


4.804 (5) 


0.76802 (3) 


0.85095 (4) 


0.006884 (2) 


0.00634 (1) 


60 


4.912 (2) 


0.86337 (4) 


0.98540 (4) 


0.014770 (3) 


0.010449 (6) 


65 


5.1(2) 


0.9723 (7) 


1.1463 (8) 


0.02653 (2) 


0.015548 (8) 


70 


5.311(6) 


1.09812 (5) 


1.34094 (6) 


0.043840 (4) 


0.02171 (1) 


75 


5.34(1) 


1.2478(2) 


1.5827(2) 


0.069408 (9) 


0.029019(10) 


80 


5.463 (2) 


1.42573(5) 


1.8834(2) 


0.107174(8) 


0.037514 (4) 


83 


5.521 (3) 


1.55007 (6) 


2.1010(2) 


0.138236 (10) 


0.043186 (7) 


90 


5.71 (5) 


1.905 (2) 


2.751 (3) 


0.2488 (2) 


0.05798 (3) 


92 


5.857(3) 


2.0238 (4) 


2.9786 (8) 


0.29400 (2) 


0.06251 (1) 


100 


5.86 (6) 


2.650 (4) 


4.224(8) 


0.5871 (8) 


0.08183 (6) 



B. Two- loop vacuum polarization 

The two-loop VP correction is represented by the dia- 
grams in Fig. [TJg)-(i). (The term the "two- loop VP" will 
be abbreviated as "VPVP" in the following.) It will be 
convenient to split our evaluation of the VPVP correc- 
tion into two parts, considering separately the diagram 
(g) and the remaining two diagrams (h) and (i). 

1. Diagram (g) 

The correction induced by the diagram in Fig.QJg) can 
be regarded as the second-order perturbation contribu- 
tion induced by the one- loop VP potential [/ypj 

a jji (a|gyp|rt) (n\U VP \a) 
ZAii/vpvp,g = 2^ • > 

The numerical evaluation of this correction is relatively 
simple. It was carried out by employing the general 
scheme developed for the VP potential and described in 
Sec. II Bl The summation over the Dirac spectrum was 
performed by the dual-kinetic-balance basis set method 

The numerical values of the energy shifts induced by 
the diagram in Fig. QJg) are presented in Table |TT] for 
the n = 1 and n — 2 states of H-like ions. The re- 
sults are expressed in terms of the function F{Za) de- 
fined by Eq. (f2"Tj) . Good agreement is found with the 



previous evaluations of this corrections [20j, |22j . Our cal- 
culation accounts for the extended nuclear charge distri- 
bution (with the Fermi model employed for the Uehling 
potential and the wave functions and the spherical shell 
model, for the Wichmann-Kroll potential). The uncer- 
tainties listed in the table include the numerical error as 
well as the estimated errors due to the models of the nu- 
clear charge distribution and due to the uncertainties of 
the values of the rms radii. The estimation of errors was 
done similarly to that for the SEVP correction. 

The higher-order part of the correction can be identi- 
fied by taking into account its Za expansion of the form 

A£Vpvp, 9 = m y—) n3 i Q 50 

+ {Za) In [{Za)- 2 ] a 61 + (Za) G(Za)\ , 

(30) 

where [II El 11 

2377 a f-m 

«50 = -7T^()|,0, (31) 
g 

fl 6i = _ 225 5 i.o> (3 2 ) 

and G(Za) is the higher-order remainder, G{0a) — ciqo- 
The remainder G(Za) inferred from our numerical re- 
sults is plotted on Fig. [3] The results for hydrogen are: 
G u {a) = -0.115, G 2s (a) = -0.059, \G 2pi/2 (a)\ < 10" 4 , 
and|G 2p3/ »| <10- 4 . 
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TABLE II: Energy shifts induced by the VPVP diagram in Fig. QJg), in units of F(Za). The uncertainties specified include 
the estimated errors due to the nuclear charge distribution models and due to the rms radii. 



z 


Is 


2s 


2pi/ 2 


2^3/2 


1 


-0.000490 


-0.000487 






2 


-0.001018 


-0.001006 






3 


-0.001578 


-0.001552 






5 


-0.002784 


-0.002716 






7 


-0.004092 


-0.003967 


-0.000001 




10 


-0.006231 


-0.005995 


-0.000002 




15 


-0.010260(1) 


-0.009803 (1) 


-0.000009 




20 


-0.014896(1) 


-0.014216(1) 


-0.000026 


-0.000001 


25 


-0.020225 (2) 


-0.019375 (2) 


-0.000060 


-0.000002 


30 


-0.026376 (3) 


-0.025480 (3) 


-0.000124 


-0.000003 


35 


-0.033524 (5) 


-0.032796 (5) 


-0.000235 


-0.000005 


40 


-0.041928 (8) 


-0.041708 (8) 


-0.000421 


-0.000007 


45 


-0.05185(1) 


-0.05265 (1) 


-0.000727 


-0.000010 


50 


-0.06374(2) 


-0.06631 (2) 


-0.001218(1) 


-0.000015 


55 


-0.07810 (3) 


-0.08356 (3) 


-0.001998(1) 


-0.000021 


60 


-0.09568 (5) 


-0.10566 (6) 


-0.003233 (2) 


-0.000029 


UJ 


—0 1 1 73 (91 


— n 1 349 


U.UUJIO ^-L 1 


_ n nnnn^Q 

u.uuuuoy 


70 


-0.1441 (1) 


-0.1713(1) 


-0.008238 (7) 


-0.000051 


75 


-0.1785 (2) 


-0.2212 (2) 


-0.01313(1) 


-0.000067 


80 


-0.2221 (3) 


-0.2880 (4) 


-0.02094 (3) 


-0.000085 


83 


-0.2542 (4) 


-0.3393 (5) 


-0.02778 (4) 


-0.000099 


90 


-0.3522 (9) 


-0.505(1) 


-0.0543(1) 


-0.000136 


92 


-0.3865 (9) 


-0.566(1) 


-0.0658 (1) 


-0.000149 


100 


-0.583 (3) 


-0.935 (4) 


-0.1482 (6) 


-0.000209 



a -i.o 
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FIG. 3: Higher-order remainder G(Za) induced by the dia- 
gram in Fig.QIg). 



2. Diagrams (h) and (i) 

We now consider the energy shift induced by the dia- 
grams in Fig.[TJh) and (i). Its leading (in Za) part can be 



obtained within the free-loop approximation, i.e., keep- 
ing the first nonvanishing term in the expansion of the 
fermion loop in terms of the binding potential. The corre- 
sponding expression was derived long ago in the classical 
paper by Kallen and Sabry 55|. Later, it was re-derived 
by a number of other techniques [H, [H?], [H[ . According 
to the Furry theorem, corrections to the Kallen-Sabry po- 
tential are suppressed by a factor of (Za) 2 , so that this 
potential is supposed to yield a dominant contribution 
even in the medium- Z region. In the present investiga- 
tion, diagrams (h) and (i) will be approached within the 
free- loop approximation only. In order to stress this fact, 
we will use the superscript "KS" in the formulas below. 

In the free-loop approximation, the correction induced 
by the diagrams in Fig. QJh) and (i) is given by the ex- 
pectation value of the Kallen-Sabry potential, 



AE™ VP!hi = {a\V KS \a) 



(33) 



For a spherically symmetric nuclear charge distribution, 
the Kallen-Sabry potential can be conveniently written 
in the form 59] 



V KS (r) = / dr' r'p(r') [L (2m\r - r'|) - L a (2m\r + r'\)\ , (34) 
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where 

^ - r*^{(-&-^<-*) ^+(s- i-i-<!) 

+ ("5?-Jj)^ rrTl ° [8t < (1 - 1) l 



and the nuclear charge density is normalized by 
J dr p(r) = 1 . 

The results of our numerical evaluation of the energy 
shift due to the Kallen-Sabry potential are listed in Ta- 
ble [Till The numerical values presente d ag ree well with 
the results of the previous studies [HI, Ufa]- Our calcu- 
lation was carried out with employing the Fermi model 
for the nuclear charge distribution. The values for the 
rms radii and their uncertainties are listed in Table UJ 
The uncertainties specified in Table IITTI for our numerical 
results include the numerical error as well as the errors 
due to the nuclear-charge distribution model and due to 
uncertainties of the rms radii. The error due to the nu- 
clear model was estimated by taking the difference of the 
results obtained with the Fermi and the uniform model. 

The Za expansion of the Kallen-Sabry contribution 
reads 



AR 



KS 

VP VP, hi 



■ m 



(Za) 2 In [(Za)- 2 ] a 6 i + (Za) 2 G(Za)} 

(36) 



where the function G(Za) is the higher-order remainder. 
The results for the first terms of the Za expansion are 



O40 
G61 



/ 45331 25 52 , 

ttH In 2 

V 39690 63 63 

~2~' 



(37) 
(38) 
(39) 



The higher-order remainder G(Za) inferred from the 
Kallen-Sabry contribution is plotted as a function of the 
nuclear charge number Z on Fig. [3J For hydrogen, the re- 
sults for the Kallen-Sabry remainder term are: G\ s (a) = 
-2.642, G 2s (a) = -3.303, G 2pi/2 (a) = -0.263, and 

G 2p3/ » = -0-073. 

Our numerical results for the higher-order remainder 
for the total VPVP correction exhibit good agreement 



52j for the 



with the analytical values obtained in Rcf. 
normalized difference of the 2s and Is states and for the 
fine-structure difference. Indeed, for the VPVP remain- 
der, our calculation yields: G2 S (a) — Gi s (a) = —0.605 



CD 




10 20 30 40 50 60 70 80 90 100 

Nuclear charge number Z 

FIG. 4: Higher-order remainder G(Za) induced by the 
Kallen-Sabry contribution. 



and G2 P3/2 (a) — G2 Pl/2 (a) — 0.190, to be compared with 
the analytical results of ago (2s) — a6o(ls) = —0.611365 
and d6o(2p3/2) — a,6o(2pi/ 2 ) = 0.189815, correspondingly. 

A complete evaluation of the VPVP correction beyond 
the free-loop approximation is a difficult problem [espe- 
cially, for diagram (i)] and has not been carried out up 
to now. For diagram (h), the calculation is easier and 
can be performed by a generalization of methods devel- 
oped for the one-loop VP correction. For uranium and 
lead, such a calculation was reported in Ref. [23]. The 
numerical values obtained for this diagram for the con- 
tribution beyond the free-loop approximation turned out 
to be rather small; it is expected that the corresponding 
contribution from diagram (i) is much larger. In the ab- 
sence of a direct calculation, we estimate the theoretical 
uncertainty of the VPVP correction due to the omitted 
part beyond the free-loop approximation by multiplying 
the absolute value of the Kallen-Sabry contribution by a 
factor of (Za) 2 . 



C. Vacuum-polarization insertion into the 
self-energy photon line 

In this section, we address the correction induced by 
the SE diagram with the VP insertion into the photon 
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TABLE III: Energy shifts due to the Kallen-Sabry potential, in units of FiZa). The uncertainties specified include the estimated 
errors due to the nuclear charge distribution models and due to the rms radii. 



z 


Is 


2s 


2pi/ 2 


2^3/2 


1 


-1.002032 


-1.002067 


-0.000014 


-0.000004 


2 


-0.992369 


-0.992509 


-0.000056 


-0.000015 


3 


-0.983258 


-0.983568 


-0.000125 


-0.000034 


5 


-0.966493 


-0.967341 


-0.000344 


-0.000093 


7 


-0.951435 


-0.953071 


-0.000670 


-0.000179 


10 


-0.931629 


-0.934898 


-0.001360 


-0.000354 


15 


-0.904993 (1) 


-0.912148 


-0.003055 


-0.000757 


20 


-0.885141 (1) 


-0.897624 (1) 


-0.005471 


-0.001287 


25 


-0.871198(1) 


-0.890494 (1) 


-0.008683 


-0.001927 


30 


-0.862607(2) 


-0.890315 (2) 


-0.012801 


-0.002669 


35 


-0.859040 (3) 


-0.896935 (3) 


-0.017978 


-0.003504 


40 


-0.860403 (4) 


-0.910519 (5) 


-0.024416 


-0.004426 


45 


-0.866633 (6) 


-0.931345 (7) 


-0.032384 


-0.005431 


50 


-0.877972 (8) 


-0.960119 (9) 


-0.042236 


-0.006517 


55 


-0.89473(1) 


-0.99776 (1) 


-0.054441 


-0.007684 


60 


-0.91747(2) 


-1.04561 (2) 


-0.069624(1) 


-0.008931 




— n Q4fS8 
u.y±uo v*-*/ 


— 1 1 D5^ (A\ 


U.UOOUi \*-) 


— n m 095Q 


70 


-0.98346 (4) 


-1.17895 (4) 


-0.112569 (2) 


-0.011673 


75 


-1.02957(7) 


-1.27069 (9) 


-0.143102 (5) 


-0.013175 


80 


-1.08601 (7) 


-1.38369 (9) 


-0.182409 (7) 


-0.014770 


83 


-1.12595 (8) 


-1.4642 (1) 


-0.21150(1) 


-0.015775 


90 


-1.2402(5) 


-1.6988 (8) 


-0.30171 (8) 


-0.018267 


92 


-1.2784(2) 


-1.7791 (2) 


-0.33488 (3) 


-0.019020 


100 


-1.476 (1) 


-2.201 (2) 


-0.5202 (3) 


-0.022227 



line, depicted by Fig. [Ijk). The corresponding shift of 
the energy will be referred to as the S(VP)E correction. 
The leading (in Za) part of this correction can be again 
obtained within the free-loop approximation. Unlike the 
VPVP contribution, however, corrections to the free-loop 
approximation in this case are suppressed only by the 
first power of Za. The leading term beyond this approx- 
imation is known from perturbative calculations |49t [fiH . 
An all-order calculation of the S(VP)E correction beyond 
the free-loop approximation is a difficult problem and has 
not been performed so far. In the present investigation, 
we will approach the S(VP)E correction within the free- 
loop approximation only. 

The evaluation of the S(VP)E correction within the 
free-loop approximation can be performed by a gener- 
alization of the method for the one-loop SE correction 
described in Section ITAl The S(VP)E correction is rep- 
resented by a sum of the zero-, one-, and many-potential 
terms: 

A£ SV PE,a = A£$ PE>a + A£« PE a + A25< a +> E a , (40) 

where the subscript "a" indicates the free-loop approxi- 
mation. Formulas for the three terms above can be ob- 



tained from the corresponding expressions in Section II Al 
by replacing the standard photon propagator by the 
"dressed" one, obtained by inserting the renormalizcd 
one-loop VP tensor into the photon line. In D = 4 — 2e 
dimensions, the replacement is given by [2ll . Ir35j 

1 _aC 1 2e f 1 z 2 {l-z 2 /i) 

k 2 +i0^ 4tt m J Z [m 2 -k 2 {l- z 2 )/A-iOY+* ' 

(41) 

where C e = (4tt) £ T(l + e) n 2e /m 2e . 

The zero-potential term is represented by 

A 4vPE,a = / j0jS**(p)ri&,B.(ea,p)1>a(p), (42) 

where the free dressed SE operator Syp is obtained from 
Eq. © by the substitution (|4"Tj) , and its renormalizcd 

part Syp R is the finite part of the difference Syp — 5m 
when e approaches zero. Applying the standard tech- 
nique for the evaluation of momentum integrals (see, e.g., 

Appendices of Ref. 28]), the difference Eyp — 8m is con- 
veniently represented as 
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,(0) 
J VP 



(p) — Sm 



47T 



2 r 



2 
i 



13 



599 
~54~ 



Stt 2 
9 



— m) 



— ] 8 / dx dz 



z 2 (l - z 2 /3) 



[(1 - x)j> - 2m] In 



x(l-z 2 )F + 4(l-x) 
x 2 (l - z 2 ) +4(1 - x) 



(43) 



where V = x(l — p) + p and p = (to 2 — p 2 )/m 2 . 
The one-potential term is given by 

A4vp E , a = y p^s^i 

X r VP , R ( £ a , Pi ; £a , P2 ) Vtf ( <? ) tpa 0»2 )(M) 



where r VP R is the finite (when e — ► 0) part of the time 
component of the dressed vertex operator r vp , which is 
obtained from the one-loop vertex operator P, Eq. (jlip . 
by the substitution (141|) . Evaluating the momentum in- 
tegrations, we obtain the following representation for the 
operator r yp , 



r VP (Pi,P2) 



47T 



2 

3? 



13 

9e 



cfa; <ij/ dz xz 



877 
"54~ 



z 



16tt' 



4iV Q (x6) 



;(1 - z 2 )L> + 4to 2 (1 - x) 1 - z s 



In 



z(l-z 2 )£>/m 2 + 4(l-a;) 
x 2 (l-z 2 )+4(l-x) 



(45) 



where I? = xb 2 +m 2 — yp\ — (1 — y)p2, b = ypi + (l — y)p2, 
and N a (k) = ^ ($t - ft + m)j a (fa - ft + m)^ a . Using the 
Ward identity, it is easy to check that the divergent parts 
of Eqs. ([IS"]) and (|4"3"j) (i.e., terms ~ e" 1 and e~ 2 ) cancel 
each other in the total expression for the energy shift. 
This is the justification of our definition of the renormal- 
ized parts of the operators £y P and Ty P , which consists 
just in dropping out the divergent part of Eqs. f4"3"|) and 

The numerical evaluation of the zero- and one- 
potential terms is similar to that for the one-loop SE 
correction but is more time consuming, due to the pres- 
ence of additional integrations over Feynman parameters. 
In order to achieve high numerical accuracy for the one- 
potential term in the \ow-Z region, we had to identify the 
contribution of Eq. (|45|) at p\ — P2 and evaluate it sep- 
arately. The subtraction of this contribution in Eq. ([44]) 
makes the integrand to be a smooth function at q = 0, 
which simplifies numerical integrations considerably. 

The many-potential term is given by 
A-Bgyp Ea =2ia / duj dxidx 2 Dy P (uj,Xi 2 ) 

J — oo J 

x ipl(x{) aip G( 2+ ) (e a - u>, x 1 , x 2 ) a v i/j a {x 2 ) 

(46) 

This formula differs from the expression for the many- 
potential part of the one-loop SE correction only by the 



dressed photon propagator D vp , which reads 



D^ p (uj,x 12 ) 



a 

7T 



dt y/t 2 - 1 



2t 2 



1 



3t 4 

x D^(Lj,x 12 ;2mt) . (47) 



D^ v {uj 1 x;X) is the propagator of the photon with a mass 
A, whose expression in the Feynman gauge is 



D*>, x 12 ; A) = fT eXp[Wuj2 ~ A2 + <0 M . (48) 

47T X\2 

It is easy to see that the numerical calculation of the 
many-potential term falls naturally into two steps: (i) the 
evaluation of the many-potential part of the one-loop SE 
correction with an effective photon mass A = 2mi and 
(ii) the numerical integration over t as given by Eq. (|47p . 
There is even a certain simplification as compared to the 
one-loop case. It is a common approach to deform the 
contour of the u> integration in Eq. (|46]) . separating it 
into the low-energy and the high-energy part (see, e.g., 
Ref. [H for details). In the case of the S(VP)E cor- 
rection, the contribution induced by the low-energy part 
of the contour vanishes identically, which is due to the 
condition on the effective photon mass A > 2m. 

The results of our numerical evaluation of the S(VP)E 
correction for the n = 1 and n — 2 states of H-like ions are 
presented in Table [TV] Our numerical results are in good 
agreement with the data obtained previously for the Is 
state [20l[2ll and with the 2s and 2pi/ 2 values for Z = 92 



from Ref. 




1 ' 2 

Our calculation was performed within the 
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free-loop approximation and for the point nuclear model. 
The uncertainty specified in the table is the numerical 
error only. We estimate the theoretical uncertainty due 
to uncalculated terms beyond the free-loop approxima- 
tion by multiplying the absolute value of the correction 
by a factor of 3 (Za). This factor arises as a ratio of 
the leading-order contribution beyond the free-loop ap- 
proximation for the Is state, —0.386 m(a/ir) 2 (Za)^ [49j |. 
and the leading-order contribution within this approxi- 
mation, 0.142 m(a/n) 2 (Za) 4 0. 

The inclusion of the finite nuclear size (FNS) effect is 
not necessary at present for the S(VP)E correction, since 
it is expected to yield a much smaller contribution than 
the error due to the free-loop approximation. The rela- 
tive contribution of the FNS effect on the S(VP)E cor- 
rection can be estimated by taking the relative values of 
this effect for the one-loop SE correction. To the leading 
orders in Za and lni?, the relative value of the FNS-SE 
effect for the ns states is [6l|, |62( 



Sfns 



- a 



4 In 2 



(49) 



where b = exp[l/(2 7 ) - C - 5/6] , 7 = y/l - (Za) 2 , C « 
0.557 is the Euler constant, and R is the nuclear radius. 
Analogous formulas for the npi/2 and np^/ 2 states can be 
found in Ref. [6^|. Numerical values of 8fns are within 
3% for the whole region of the nuclear charge numbers. 

The Za expansion of the S(VP)E correction within the 
free-loop approximation reads 



AE. 



SVPE ,a 



(Za) 4 



|a 40 + (Za) a 50 



+ (Za) 2 In [(Za)- 2 ] a el + (Za) 2 G(Za)} , 

(50) 

where the function G(Za) is the higher-order remainder. 
The results for the first terms of the Za expansion are 



O40 



7 

81 



119 
"36" 



5tt 2 
216 

_7T 2 

3 



3(3 + 1) - 1(1 + 1) - 3/4 



I(I+1)(2J+1) 



(1 



a 50 = -0.22 9 5 3 6 W , 

a 40 r 

~2~ 



061 



J 10 ■ 



8l0), 

(51) 
(52) 

(53) 



The higher-order remainder G(Za) inferred from our 
numerical results is plotted in Fig. [5l For hydrogen, 
our results are: G ls (a) = 0.93(6), G 2s (a) = 1.04(5), 
G 2pi/2 (a) = 0.02(4), and G 2p3/2 (a) = 0.01(3). These 
values are consistent with the Za-expansion results for 
the normalized difference of the 2s and Is states and 



TABLE IV: Energy shifts due the S(VP)E correction eval- 
uated within the free-loop approximation and for the point 
nucleus, in units of F(Za). 



Is 



2 s 



2pi/ 2 



10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

60 

65 

70 

75 

80 

83 

90 

92 

100 



140459 (4) 

139002 (4) 

137650 (4) 

135219 (4) 

133099 (3) 

13042 (1) 

12708(1) 

12493(1) 

12381 

12366 

12444 

12615 

12882 

13254 

13742 

14362 

15137 

16099(1) 

17291 (1) 

18773(1) 

19835 (1) 

22987(2) 

24110(2) 

30064 (3) 



0.140465 (3) 

0.139026 (2) 

0.137705 (2) 

0.135370(2) 

0.133390 (4) 

0.13100(1) 

0.12835(1) 

0.12715(1) 

0.12727 

0.12866 

0.13134 

0.13536(1) 

0.14085(1) 

0.14801(1) 

0.15711(1) 

0.16851 (1) 

0.18271 (1) 

0.20042 (1) 

0.22260(1) 

0.25059(2) 

0.27097(2) 

0.33285 (3) 

0.35535 (3) 

0.47842 (6) 



-0.005228(2) 

-0.005226(1) 

-0.005223(1) 

-0.005212(1) 

-0.005196(3) 

-0.00516 

-0.00507 

-0.00492 

-0.00471 

-0.00442 

-0.00403 

-0.00350 

-0.00280 

-0.00187 

-0.00064 

0.00098 

0.00314 

0.00603 

0.00994 

0.01528 

0.01943 

0.03322 

0.03863 

0.07124(1) 



2^3/2 



0.002615(2) 

0.002615(1) 

0.002615(1) 

0.002614(1) 

0.002612 (3) 

0.00260 

0.00259 

0.00258 

0.00255 

0.00252 

0.00249 

0.00246 

0.00242 

0.00237 

0.00232 

0.00226 

0.00220 

0.00213 

0.00205 

0.00197 

0.00191 

0.00177 

0.00172 

0.00152 
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FIG. 5: Higher-order remainder G(Za) for the S(VP)E cor- 
rection within the free-loop approximation. 



for the fine-structure difference [52j: a 60 (2s) — a 6 o(ls) = 
0.109999 and a 60 (2p 3/2 ) - a 60 (2p 1/2 ) = -0.013435. 

For completeness, we specify also the result for the 
leading term of the Za expansion known for the S (VP) E 
correction beyond the free- loop approximation [49L l64| , 



A£ SV PE,b = m [-0.38615 6 lQ ] 



(54) 
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III. CONCLUSION 

In the present investigation, we performed calculations 
of the part of the two-loop Lamb shift induced by the di- 
agrams in Fig.[TJd)-(k). Numerical results were obtained 
for the n = 1 and n = 2 states and for the whole region 
of the nuclear charge numbers Z = 1 — 100. The dia- 
grams (d)-(g) were calculated rigorously to all orders in 
Za, whereas for the diagrams (h)-(k), the fermion loops 
were approximated by their leading Za-expansion con- 
tribution. An estimate was given for the higher-order 
terms thus omitted. The finite nuclear size effect was ac- 
counted for in the evaluations of all diagrams except for 
the diagram (k). The latter diagram was calculated with 
the point nuclear model; an estimate of the finite nuclear 
size effect was supplied. 

In the low-Z region, our numerical results were em- 
ployed for the identification of the nonperturbative re- 
mainder G(Za), which incorporates all orders in the 
Za-expansion starting with a 2 (Za) e . For hydrogen, the 
net result for the two-loop diagrams with closed fermion 
loops is 



G ls (a) = 
G 2 » = 
G 2pi/2 (a) = 



G 



2p 3 /2 



(a) 



-15.0(4) (2.2), 
-14.0(4) (2.2), 
-0.28(4) , 
-0.05(3), 



(55) 
(56) 
(57) 
(58) 



where the first error quoted is the numerical uncertainty. 
The second error (if given) is due to contributions beyond 
the free-loop approximation. We estimate them for the 
ns states as 1 [in units of G(Za)] for the VPVP correc- 
tion and as 2 for the S(VP)E diagram (which arises as a 
typical coefficient of 0.2 enhanced by the first power of 



logarithm). It is interesting to note that the dominant 
part of the remainder term for the ns states is due to the 
SEVP correction [Fig. H(cl)-(f)]. 

In order to get the complete results for the two-loop 
Lamb shift, one should combine the numerical values ob- 
tained in the present work with the contribution due to 
the two-loop self-energy [Fig. [lja)-(c)]. Its all-order cal- 
culation was accomplished in our previous investigations, 

in Refs.f2im,[29|,[33] for the ls state and z > 10 and 
in Ref. [31| for the n = 2 states and Z > 60. Combined 
together, our calculations yield results for the total two- 
loop QED correction, which improve the total theoret- 
ical accuracy of the ls Lamb shift [27| and that of the 
2 Pi/2,3/2 - 2s transition energy in heavy Li-like ions 

Still, the project of the calculation of the two-loop 
Lamb shift is far from being finished. There are sev- 
eral reasons for this. First, the results of the all-order 
calculation of the two-loop self-energy correction for the 
ls state in the low-Z region |29| do not agree well with 
the Za-expansion result of Ref. [lj| . Second, the calcula- 
tion [3l[ for the n = 2 states was performed in the high- 
Z region only. Third, a part of the two-loop diagrams 
with closed fermion loops is presently calculated within 
the free-loop approximation only. Each of these points 
represents a difficult problem and all of them should be 
solved before the calculation of the Lamb shift to order 
a 2 is completed. 
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